#ifndef cardinal_spline_h
#define cardinal_spline_h
#include "spline.h"
using namespace std;
class cardinal_spline{
private:
    Function &Func;
    Point X;
    int order;
    MatrixXf solve_axb(MatrixXf A,MatrixXf b){
        MatrixXf x;
        x=A.fullPivLu().solve(b);
        return x;
    }
public:
    cardinal_spline(Function &Func,Point X):Func(Func),X(X){};
    void condition(string z,Point X,int n,double xx){
        Point Z;
        if (z=="complete_cardinal_quadratic"){
            //already correct
            MatrixXf b(n,1);
            MatrixXf A(n,n);
            MatrixXf a(n,1);
            for(int j=0;j<n;j++){
                if(j==0){
                    b(j,0)=8*Func(X.x[j])-2*Func(1);
                }
                else if (j==n-1){
                    b(j,0)=8*Func(X.x[j])-2*Func(n+1);
                }
                else b(j,0)=8*Func(X.x[j]);
            }
            for(int j=0;j<n;j++){
                if(j==0){
                    A(j,0)=5;
                    A(j,1)=1;
                }
                else if(j==n-1){
                    A(j,j-1)=1;
                    A(j,j)=5;
                }
                
                if(j!=0&&j!=n-1){
                    A(j,j)=6;
                    A(j,j-1)=1;
                    A(j,j+1)=1;
                }
            }
            a=solve_axb(A, b);
            /*cout<<A<<endl;
            cout<<b<<endl;
            cout<<endl<<a<<endl;*/
            double result=0;
            if(xx<0){
                cout<<xx;
                xx=-xx;
            }
            else{
                cout<<xx;
            }
            for(int i=0;i<=n+1;i++){
                if(i==0){
                    result+=(2*Func(1)-a(0,0))*Bspline(Func, X, 2).B_i(xx, i, 2);
                }
                else if(i==n+1){
                    result+=(2*Func(n+1)-a(n-1,0))*Bspline(Func, X, 2).B_i(xx, i, 2);
                }
                else {
                    result+=a(i-1,0)*Bspline(Func, X, 2).B_i(xx, i, 2);
                }
            }
            cout<<" "<<result<<endl;
        }
        else if(z=="complete_cardinal_cubic"){
            //already correct
            MatrixXf b(n,1);
            MatrixXf A(n,n);
            MatrixXf a(n,1);
            for(int j=0;j<n;j++){
                if(j==0){
                    b(j,0)=3*Func(X.x[j])+Func.diff(X.x[j]);
                }
                else if (j==n-1){
                    b(j,0)=3*Func(X.x[j])-Func.diff(X.x[j]);
                }
                else b(j,0)=6*Func(X.x[j]);
            }
            for(int j=0;j<n;j++){
                if(j==0){
                    A(j,0)=2;
                    A(j,1)=1;
                }
                else if(j==n-1){
                    A(j,j-1)=1;
                    A(j,j)=2;
                }
                if(j!=0&&j!=n-1){
                    A(j,j)=4;
                    A(j,j-1)=1;
                    A(j,j+1)=1;
                }
            }
            a=solve_axb(A, b);
            /*cout<<endl<<A<<endl;
            cout<<b<<endl;
            cout<<endl<<a<<endl;*/
            if(xx<0){
                cout<<xx;
                xx=-xx;
            }
            else cout<<xx;
            double result=(a(1,0)-2*Func.diff(X.x[0]))*Bspline(Func, X, 3).B_i(xx, -1, 3);
            for(int i=0;i<=n;i++){
                if(i==n){
                    result+=(2*Func.diff(X.x[n-1])+a(n-2,0))*Bspline(Func, X, 3).B_i(xx, i, 3);
                }
                else {
                    result+=a(i,0)*Bspline(Func, X, 3).B_i(xx, i, 3);
                }
            }
            cout<<" "<<result<<endl;
        }
    }
};
#endif /* cardinal_spline_h */
